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1.  INTRODUCTION 


This  final  report  consists  of  a  summation  of  results  of  a 
survey  conducted  by  Control  Dynamics  Co.  from  December  1984 
through  February  1985.  The  purpose  of  this  survey  was  to 
determine  the  type  of  work  being  done  in  the  field  of  wind 
and  turbulence  modeling  as  related  to  the  problems 
encountered  In  the  use  of  free-fllght  flat-fire  rockets  In 
the  atmospheric  boundary  layer.  In  addition  it  was  desired 
that  the  modeling  techniques  which  seemed  particularly 
relevant,  be  explained  from  an  applied  engineering  point  of 
view. 

The  first  section  of  this  report  contains  a  summary  of 
surveys  of  the  available  literature  relevant  to  modeling  low 
level  winds.  The  second  section  to  this  report  consists  of 
a  summary  of  meetings  and  phone  contacts  between  Control 
Dynamics  and  various  authorities  In  the  fields  of  atmos¬ 
pheric  science.  Section  3  consists  of  definitions  of  some 
commonly  encountered  terms  used  to  described  the  assumptions 
used  in  modeling  turbulence.  Sections  5  through  8  consist 
of  discussions  of  specific  work  which  appears  relevant  to 
the  modeling  of  wind  with  regard  to  rocketry  problems. 

These  topics  Include  the  approximation  of  the  von  Karman 
turbulence  spectra  with  a  rational  function,  effects  of 
atmospheric  stability  on  turbulence,  data  available  on  mean 
winds  and  standard  deviations,  and  the  spatial  modeling  of 
wind  and  turbulence. 


2.  SURVEY  OF  AVAILABLE  LITERATURE 

A  survey  of  literature  relevant  to  boundary  layer  wind  effects 
on  flat  fire  rockets  was  conducted  through  several  channels. 
The  majority  of  the  survey  consisted  of  literature  searches 
through  the  Redstone  Scientific  Information  Center  (RSIC)  of 
the  following  data  bases. 

1.  Defense  Technical  Information  Center  (DTIC) 

2 .  NASA  Recon 

3.  Journals  and  reference  material 

Additional  information  was  obtained  from  direct  contacts  with 
several  people  working  in  the  field  of  wind  and  turbulence 
research. 

Approximately  80-90  papers  have  been  identified  as  having 
some  degree  of  relevance  to  this  task.  These  documents  fall 
roughly  into  the  following  categories. 

1.  Handbooks 

2.  Terrain  effects 

3.  Space  -  time  relationships 

4.  Flat  -  fire  trajectories 

5.  Wind  profiles,  studies,  ect. 

6.  Modeling 

7.  General  turbulence  analysis 

8.  Mi sc. 

The  majority  of  these  deal  with  the  simulation,  modeling,  or 
representation  of  various  types  of  low  -  level  turbulence. 

All  documents  obtained  by  Control  Dynamics  have  been  listed 
In  an  extensive  bibliography  In  the  last  part  of  this  report. 
Because  of  the  number  of  reports,  the  bibliography  has  been 
broken  down  into  sections  roughly  according  to  topic  as  shown 
above.  All  references  in  this  report  will  have  the  following 
form. 

(category  number  .  document  number) 

For  example,  (1.1)  would  refer  to  the  document  in  the 


3.  MEETINGS  AND  PHONE  CONTACTS 


Several  meetings  and  phone  conversations  have  been  held 
regarding  the  wind  project.  These  are  listed  in  roughly 
chronological  order  with  the  highlights  of  each  contact. 


MEETINGS 

1.  Dr.  Oscar  Essenwanger  12-4-84  MICOM  -  Research 
Di rectorate 

Discussed  the  availability  of  instrumentation  capable 
of  measuring  at  high  frequencies  (>  1  Hz.).  Best  data 
available  from  his  office  is  at  1  second  but  is  all 
very  low  speed  (<  4  m/s).  Next  best  data  is  at  6 
second  intervals.  Dr.  Essenwanger  considers  terrain  to 
be  the  most  important  factor  in  low  level  turbulence. 

2.  Don  Combs  and  Dick  Dickson  -  12/5/84  -  MICOM 

Discussed  objectives  of  the  contract  and  the  delivered 
scope  of  work.  Basic  objective  was  for  a  general  study 
which  did  not  focus  on  a  particular  missile  system.  It 
was  desired  to  look  at  the  flat  fire  scenario  in  the 
atmosphere  below  200  meters.  Material  was  reviewed  and 
received  from  D.  Dickson  on  some  basic  sources  of  infor¬ 
mation. 

3.  Dr.  Warren  Campbell  1-8-85  Systems  Dynamics  Lab  (MSFC/NASA) 

A  meeting  was  held  with  Dr.  Campbell  to  discuss  several 
issues  with  which  he  is  knowledgeable.  These  include  the 
Space  Shuttle  Turbulence  Tapes  (SSTT),  rational  approxi¬ 
mation  of  the  von  Karman  turbulence  spectra,  and  the 
inclusion  of  cross  correlation  factors  into  turbulence 
models.  The  SSTT  are  available  for  use.  Several  points 
regarding  their  applicability  were  discussed  as  follows: 

1.  The  SSTT  are  generated  over  a  frequency  range 
based  on  the  size  of  the  vehicle  in  question 
(Space  Shuttle),  thus  to  be  applied  to  short 
tactical  rockets  the  tapes  would  need  to  be 
regenerated.  If  it  was  desired  to  use  the  tapes 
with  larger  strategic  missiles,  the  present 
tapes  would  seem  applicable. 

2.  SSTT  data  exists  im  six  separate  altitude  bands 
and  was  generated  with  differing  sample  rates  in 


various  bands.  If  the  tapes  were  to  be  used  for 
simulation  of  lofted  missile  trajectories  which 
pass  through  more  than  one  band,  the  variations 
in  sample  rate  would  have  to  be  allowed  for. 

3.  Use  of  the  tapes  in  simulation  work  requires  the 
use  of  a  computer  with  sufficient  storage  capa¬ 
bility  to  store  the  entire  turbulence  time  series 
needed  for  a  simulation  run  in  a  file.  The  simu¬ 
lation  would  then  read  data  from  mass  storage  as 
it  ran. 

4.  The  SSTT  time  series  have  been  analyzed  statisti¬ 
cally  and  have  been  shown  to  possess  very  good 
spectral  characteristics.  The  tapes  exhibit  a  -5/3 
rolloff  as  associated  with  the  von  Karman  spectra. 

A  paper  was  received  from  Dr.  Campbell  regarding  a 
method  for  approximating  the  von  Karman  spectra  by  a 
rational  expression.  The  technique  used  shows 
extremely  good  approximation  to  the  theoretical  von 
Karman  spectra.  Use  of  an  approximation  would  allow 
the  modeling  of  the  more  accurate  von  Karman  spectra 
without  the  associated  difficulty  of  simulating  the 
irrational  von  Karman  expressions.  Dr.  Campbell 
delivered  two  computer  programs  to  Control  Dynamics 
which  use  the  approximate  von  Karman  techniques  to 
produce  turbulence  time  series.  An  elaboration  on  this 
work  will  be  presented  in  a  following  section. 

A  technique  which  is  currently  under  study  by  Dr. 

Campbell  was  presented  for  the  inclusion  of  cross 
correlation  terms  into  the  turbulence  models.  The 
technique  involves  the  process  of  filtering  various 
white  noise  inputs  and  summing  the  outputs  through 
filters  to  produce  turbulence  with  the  desired  cross 
correlation  and  spectra. 

Two  additional  sources  of  information  were  suggested 
as  follows: 

1.  Models  have  been  developed  by  Paul  Reeves  to 
simulate  the  patchy  or  modulated  characteristics 
of  turbulence. 

2.  It  was  suggested  that  Dr.  Walter  Frost  in  Tullahoma 
Tenn.  (UTS I  -  University  of  Tennessee  Space  Insti¬ 
tute)  be  contacted  for  a  knowledgeable  source  of 
information  on  both  flight  mechanics  and  turbulence 
model i ng . 
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4.  Dr.  Oscar  Essenwanger  2-11-85  MICOM  -  Research 
Directorate 

A  follow-up  meeting  was  held  with  Dr.  Essenwanger  to 
clarify  information  received  in  a  meeting  on  12-4-84. 

It  was  also  desired  to  obtain  specific  information  on 
the  following  topics. 

a.  Mean  wind  profiles  and  availability  of  data  for 
the  lower  atmosphere. 

b.  Information  on  relative  stability  of  the 
atmosphere. 

c.  Availability  of  any  turbulence  or  any  other  wind 
model . 

d.  Information  on  any  work  done  on  correlation  and 
coherence  of  turbulence. 

e.  Any  information  relating  to  terrain  Influences 
on  wind  and/or  turbulence. 

f.  Any  information  concerning  the  non-gaussian 
nature  of  turbulence. 

Work  done  by  Dr.  Essenwanger* s  office  has  been  primarily 
limited  to  climatological  studies  with  very  little  work 
done  on  turbulence.  Although  two  reports  have  been 
received  concerning  turbulence  representation,  no  models 
are  available  from  Dr.  Essenwanger.  The  reports  that 
are  available  do  not  seem  to  be  useful  In  the  current 
modeling  effort. 

No  significant  work  has  been  performed  by  Dr.  Essenwanger* s 
office  in  the  area  of  terrain  effects,  correlation,  coher¬ 
ence,  or  the  non-gaussian  nature  of  turbulence. 

Two  areas  of  work  from  Dr.  Essenwanger* s  office  appear 
to  have  relevance  to  low-level  wind  modeling.  The  first 
area  consists  of  description  of  mean  winds  encountered 
for  various  seasons  and  geographic  locations.  This  data 
is  available  from  Dr.  Essenwanger  in  the  form  of  numerical 
tables.  Work  has  also  been  published  compiling  studies  of 
atmospheric  stability  observed  in  Frankfurt  Germany  and 
Osan  Korea.  This  study  should  provide  Insight  Into  the 
relative  occurrence  of  stable,  unstable,  and  neutral  condi¬ 
tions  in  the  atmosphere.  This  information  will  aid  under¬ 
standing  the  area  of  applicability  of  various  turbulence 
models. 
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Dr.  Warren  Campbell  2-22-85  Systems  Dynamics  Lab 
(MSFC/NASA) 

A  second  meeting  was  held  with  Dr.  Campbell  to  discuss 
the  spatial  nature  of  turbulence.  Various  approaches 
for  Incorporating  coherence  Into  turbulence  models  were 
discussed.  Dr.  Campbell  is  currently  working  on 
methods  to  model  turbulence  with  coherence  by  the  use 
of  Monte  Carlo  techniques  using  linear  filters.  These 
techniques  show  promise  of  providing  turbulence 
simulations  that  approximate  the  von  Karman  spectra 
and  Include  coherence. 

Also  discussed  were  the  relative  merits  of  Taylor's 
hypothesis  to  describe  the  spatial  characteristics  of 
the  wind.  It  has  been  shown  that  over  uniform  terrain 
Taylor's  hypothesis  Is  valid. 

Dr.  Walter  Frost  2-25-85  UTSI  -  FWG  Associates 

A  meeting  was  held  with  Dr.  Frost  In  order  to  discuss 
his  work  in  the  area  of  turbulence  and  atmospheric 
science  research.  Dr.  Frost  has  had  a  wide  range  of 
experience  In  the  field  of  atmospheric  science  and 
aerodynamic  research.  This  work  Includes  several 
programs  In  which  experimental  data  was  taken  either 
with  tower  arrays  on  the  ground  or  using  Instrumented 
aircraft.  Dr.  Frost  also  has  developed  extensive 
simulation  capabilities.  The  following  list  contains 
a  brief  summary  of  projects  and  capabilities  available 
from  Dr.  Frost. 

1.  Linear  filter  turbulence  simulations. 

2.  Non-gausslan  turbulence  simulation. 

3.  Inter-level  coherence  modeling. 

4.  Diffusion  models  (particle  tracking). 

5.  Models  based  on  Kalmal  spectra. 

6.  Test  data  taken  at  MSFC  tower  array  to  test 
the  validity  of  Taylor's  hypothesis. 

7.  Work  simulating  wind  shear  and  microbursts. 

8.  Data  was  collected  at  the  MSFC  tower  array. 


9.  Work  has  been  performed  with  the  NASA  gust 
gradient  program  using  an  instrumented  B-57 
aircraft. 

10.  Studies  of  complex  terrain  using  water  tunnel 
simulation  and  video  taping  facilities. 

11.  LIDAR  turbulence  measurements  in  conjunction 
with  the  B-57  aircraft  program. 

12.  Turbulence  probability  distribution  studies  to 
determine  non-gausslan  nature  of  atmosphere. 

13.  Flight  simulator  modeling  of  turbulence  and  mean 
wind.  Turbulence  and  mean  wind  models  for 
various  types  of  simulation. 

Dr.  Frost  has  accumulated  experience  in  simulating 
various  aspects  of  turbulence  such  as  occurrence  of 
various  spectra,  coherence,  terrain  effects,  ect.  Use 
of  Dr.  Frost  in  a  consulting  role  for  future  work  in 
the  area  of  developing  more  accurate  simulation 
capabilities  was  discussed. 

Dr.  Walter  Frost  2-27-85  UTSI  -  FWG  Associates 

A  follow-up  meeting  was  held  with  Dr.  Frost  at  his 
facility  In  Tullahoma  Tennessee  to  discuss  further  his 
capabilities  with  regard  to  tactical  and  other  missile 
systems.  The  simulation  and  modeling  capabilities  that 
have  been  developed  primarily  for  aircraft  systems  lend 
themselves  very  well  to  the  simulation  of  wind  effects 
on  missiles.  The  capability  to  provide  various  complex 
computer  models  for  turbulence  and  other  wind  effects 
Is  augmented  by  water  tunnel  simulation  and  video  tape 
techniques.  The  ability  to  simulate  in  the  water  tunnel 
real  terrain  models  along  with  missile  firings  and 
record  the  flow  patterns  on  video  tape  provides  insight 
Into  flow  problems  over  complex  terrains.  This  is  an 
area  which  is  very  difficult  to  address  with  numerical 
modeling  techniques.  Dr.  Frost's  water  tunnel  work  has 
been  used  in  areas  such  as  the  following: 

1.  Prediction  of  the  dispersion  of  the  Space  Shuttle 
exhaust  plume  and  ground  cloud  at  yandenburg  Air 
Force  Base.  The  water  tunnel  tape  clearly  shows 
how  the  cloud  will  dispurse  over  the  mountainous 
areas  around  the  base  for  varying  wind  speeds 
and  atmospheric  stabilities. 
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2.  Simulation  of  the  conditions  surrounding  the 
crash  of  the  Texasgulf  jet  star  aircraft  in  a 
wind  shear  related  accident.  The  water  tunnel 
correctly  predicted  the  complex  wind  flows  and 
turbulence  in  effect  at  the  time  of  the  crash 
and  was  instrumental  evidence  in  the  ensuing 
law  suits. 

3.  Simulation  of  the  exhaust  plume  released  by  a 
missile  between  launch  and  Impact  to  determine 
whether  an  obscurity  of  the  target  Is  likely  to 
occur  and  to  what  extent. 

The  combination  of  the  numerical  modeling  techniques  with 
water  tunnel  simulation  of  the  terrains  of  Interest  would 
seem  to  extend  the  insight  into  effects  of  wind  on  missile 
performance  over  that  which  Is  gained  by  numerical  models 
alone. 


PHONE  CONTACTS 

1.  Dr.  Jim  Hal  leek  DOT-Trans.  Systems  Center  617/494-2199 
Dr.  Hal leek  suggested  the  following  two  sources  of  infor¬ 
mation  on  low  level  wind  and  turbulence. 

1.  Lo-Cat  data  generated  by  the  Air  Force  during  low- 
level  wind  studies. 

2.  Boeing  document  produced  for  the  FAA. 

Barr.N. ,Gangaas,D.  and  D. Schaeffer, 1974,  Wind  Models 
for  Flight  Simulator  Certification  of  Landing  and 
Approach  Guidance  and  Control  Systems.  FAA-RD-74-206. 

2.  Dr.  A1  Bedard  NOAA-Boul der  303/497-3000  Dr.  Bedard 
Is  involved  with  the  study  of  low-level  organized  wave 
activity.  This  activity  can  produce  large  amplitude 
wavetralns  which  cause  lift  forces  on  vehicles  flying 
turough  the  wavetralns.  A  study  is  currently  concluding 
at  the  Boulder  tower  In  which  data  (taken  at  10  Hz.)  was 
collected  to  study  wave  activity.  Preliminary  information 
has  been  received  from  Dr.  Bedard  with  additional  infor¬ 
mation  to  be  available  In  approximately  the  March  1985 
time  frame. 

3.  Dr.  Walter  Frost  UTSI  615/455-0631  Dr.  Frost  is  the 
Director  of  Atmospheric  Sciences  at  the  University  of 
Tennessee  Space  Institute  in  Tullahoma  Tennessee.  Dr. 
Frost's  primary  work  has  been  in  the  area  of  wind  and 


turbulence  simulation  for  use  In  flight  simulators  and 
to  study  the  effects  of  wind  and  turbulence  on  airfoils. 
Work  has  also  been  done  In  the  area  of  terrain  effects 
and  terrain  generated  turbulence  modeling.  Much  of  the 
terrain  effects  work  Is  performed  by  using  a  water 
tunnel  to  study  flow  patterns  around  specific  terrain 
models.  Dr.  Frost  has  also  been  Invol /ed  with  the  B-57 
aircraft  turbulence  tests  and  with  the  remote  sensing  of 
turbulence  using  LIDAR. 

4.  Bill  Rodgers  -  NWS  -  772-9876  (now  retired) 

Suggested  Kelly  Hill  or  Mike  Susko  as  possible  contacts 
at  Atmospheric  Sciences  Lab  at  NASA. 

5.  Kelly  Hill  -  At.  Sc.  NASA  -  453-4175 

Suggested  Dr.  George  Flchtl  (453-0875)  or  Dr.  Warren 
Campbell  (453-1886)  as  contacts  In  the  area  of  low  - 
level  turbulence. 

6.  Dr.  Warren  Campbell  -  NASA  -  453-1886 

Dr.  Campbell  Is  primarily  concerned  with  the  spatial 
modeling  of  turbulence  In  connection  with  space  shuttle 
flight  simulation.  He  suggested  several  sources  of 
material  and  sent  several  reports  concerning  turbulence 
modeling  and  the  Space  Shuttle  Turbulence  Tapes. 

7.  FAA/ATL  -  Suggested  David  Redhun  as  a  point  of  contact. 

8.  David  Redhun  -  FAA/DCA  -  (202)426-8714 

Mr.  Redhun  Is  concerned  with  hardware  Installations  of 
low  -  level  wind  shear  measurement  equipment.  He 
suggested  three  possible  contacts  with  DOT  and  NOAA. 

-  DOT  -  Transportation  Systems  Center 
Dr.  Jim  Halleck  -  (617)494-2199 
Dr.  David  Burnham  -  (617)494-2579 

-  NOAA  -  Boulder  Tower 

Dr.  A1  Bedard  -  (303)497-6508 

9.  Jim  Bllbro  -  NASA  -  453-1597 

Mr.  Bllbro  suggested  we  contact  William  Cliff  at 
Battel le. 

10.  William  Cliff  -  Battelle  -  (509)375-2024 

Dr.  Cliff  sent  two  documents  on  the  simulation  of 
turbulence  (simple  model)  and  the  simulation  of  the 
hourly  wind  at  randomly  dispersed  sites. 


Orville  E.  Smith  -  NASA  -  453-3101 
Mr.  Smith  was  retiring  in  a  few  days  and  no  information 
was  available  from  him  because  all  of  his  references 
were  packed  away. 

Tom  Preis  (505)  678-5421 

Mr.  Preis  was  unavailable  for  comment  on  repeated  phone 
contacts. 
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A.  TURBULENCE  NOMENCLATURE 

In  any  discussion  of  atmospheric  turbulence  there  are  several 
terms  which  will  be  encountered.  These  terms  are  discussed  in 
order  to  provide  insight  into  the  information  in  following 
sections. 

1.  Homogeneity  -  A  turbulent  flow  Is  homogeneous  if  its 
statistics  do  not  vary  in  space.  The  presence  of  the 
earth's  surface  is  evidently  Important  in  two  ways  in 
considering  homogeneity.  First,  statistics  will  vary 
with  distance  from  the  boundary  so  that  It  is  unlikely 
that  homogeneity  could  prevail,  even  approximately, 
except  in  the  horizontal.  Second,  If  the  terrain  Is 
Inhomogeneous,  with  hills  and  valleys,  or  with  cities, 
fields,  and  forests,  then  the  flow  near  the  boundary 
can  hardly  be  expected  to  be  horizontally  homogeneous 
because  of  the  effects  of  the  boundary  on  the  flow 
itself.  Homogeneous  flow  can  be  valid  over  rough 
terrain  as  well  as  over  smooth  terrain  as  long  as  the 
terrain  is  consistently  rough  along  the  areas  of 
Interest.  It  Is  often  the  case  over  flat  and  homo¬ 
geneous  terrain  or  over  the  ocean  that  small-scale 
flows  seem  to  be  homogeneous(7.5). 

The  assumption  of  vertical  homogeneity  is  almost  never 
valid  near  the  ground,  the  mean  wind  and  temperature 
vary  rapidly  with  height  at  first,  and  then  somewhat 
less  rapidly  at  higher  altitudes. 

2.  Statlonarity  -  Stationarlty  Implies  that  the  statistical 
properties  of  variables  do  not  change  with  time.  For 
events  with  time  scales  of  hours  this  presents  a  problem 
due  to  the  daily  changes  of  the  atmosphere.  For  scenarios 
such  as  missile  launches  with  flight  times  of  minutes  or 
even  seconds,  the  assumption  of  statlonarity  would  seem  to 
be  quite  valid{7.5). 

3.  Isotropy  -  Isotropy  Implies  that  the  statistics  of  the 
motions  are  invariant  to  changes  in  directions  of  the 
coordinates,  and  is  usually  considered  only  when  homo¬ 
geneity  Is  present  In  all  directions.  If  Isotropy  pre¬ 
vails,  then  the  variance  of  the  three  velocity  components 
must  be  equal,  since  the  variance  of  one  component  changes 
Into  that  of  another  as  the  coordinate  system  rotates  by 
90  degrees.  But  In  fact,  the  variances  of  the  velocity 
components  In  the  boundary  layer  are  not  the  same,  and  for 
this  reason  alone,  it  Is  clear  that  the  small-scale  motions 
in  the  lower  atmosphere  are  not  Isotropic. 
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Even  though  Isotropy  does  not  apply  near  the  ground, 
local  isotropy  does.  Local  isotropy  implies  that  not 
all  the  small-scale  motions,  but  only  the  fluctuations 
on  the  very  smallest- scales,  are  isotropic.  From  this 
hypothesis  certain  statistical  properties  of  the  high- 
frequency  fluctuations  can  be  derived.  Measurements 
of  boundary-layer  turbulence  have  shown  that  these 
predictions  are  Indeed  satisfied,  provided  that  the 
size  of  the  eddies  involved  is  small  compared  to  the 
distance  to  the  surface(7.5). 
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5.  SPACE  SHUTTLE  TURBULENCE  TAPES 


One  area  which  appears  promising  is  the  use  of  the  Space 
Shuttle  Turbulence  Tapes  (SSTT)  developed  by  NASA  for  study¬ 
ing  the  effect  of  turbulence  and  gusts  on  the  space  shuttle. 
SSTT  consists  of  six  magnetic  tapes  containing  time  series 
of  various  types  of  gusts  and  gust  gradients  for  six  altitude 
bands.  The  tapes  contain  gusts  ul,  u2,  and  u3  for  the  X,  Y,  and 
Z  axes  (1,2,3)  and  gust  gradients  du2/dxl(yaw) ,  du3/dxl( pitch) , 
and  du3/dx2 ( roll ) .  The  tapes  were  produced  by  a  non-recursive 
turbulence  model  based  on  von  Karman  spectra  with  finite  upper 
limits  corresponding  to  the  dimension's  of  the  space  shuttle. 
These  upper  integration  limits  increase  with  altitude  because 
they  are  based  on  turbulence  scale  lengths  (which  Increase  with 
altitude)  and  on  characteristic  vehicle  length  in  each  axis. 

The  limits  are  calculated  according  to 

Wijmax  *  1.339  Li  /  Ij 

Li  *  turbulence  scale  length 
lj  ■  characteristic  vehicle  length 

In  typical  use  the  tapes  are  scaled  in  a  manner  considered 
appropriate  by  the  user  for  various  standard  deviation, 
intensities,  and  times.  The  dimensional  gust  is  obtained  by 
scaling  the  dimensionless  gust  by  the  appropriate  standard 
deviation  as  follows. 

ul*  =  si  ui 

ui  *  dimensionless  gust 
ui*  =  dimensional  gust 
si  *  standard  deviation 

The  dimensional  gust  gradient  is  obtained  by  scaling  the 
dimensionless  gust  gradient  by  the  appropriate  standard 
deviation  and  turbulence  scale  length  as  follows. 

dui*  si  dui 


dxj*  Li  dxj 

si  *  standard  deviation 
Li  =  turbulence  scale 

The  dimensional  time  step  is  obtained  by  scaling  the 
dimensionless  time  step  by  the  turbulence  scale  length  and 
the  vehicle  velocity  as  follows. 


LI  =  turbulence  scale 

T1  *  dimensionless  time  step 

V  =  vehicle  velocity 


Spectral  analysis  of  each  tape  reveals  that  the  simulated 
turbulence  possesses  the  appropriate  von  Karman  spectral 
characteristics.  Both  the  simulated  gusts  and  gust  gradients 
are  normally  distributed  with  near  zero  means.  The  standard 
deviation  of  each  series  is  constant  with  the  theoretical 
energy  content. (6. 28429) 

The  tapes  containing  the  turbulence  data  as  generated  for 
the  space  shuttle  are  currently  available  from  NASA.  The 
source  code  that  would  be  needed  to  regenerate  the  data  is 
not  currently  available.  New  source  code  could  easily  be 
written  by  Control  Dynamics  using  Its  FFT  (Fast  Fourier 
Transform)  subprogram. 


6.  VON  KARMAN  APPROXIMATION 


Experimental  data  has  shown  that  the  velocity  spectra  of 
atmospheric  turbulence  most  closely  corresponds  to  the 
von  Karman  spectra  which  is  of  the  following  form. 
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Figure  1.  Expressions  for  longitudinal  and  lateral  von  Karman  spectra 


It  can  be  seen  that  these  spectra  are  irrational  in  form  and 
produce  a  -5/3  rolloff  at  high  frequencies. 

Attempts  to  model  these  forms  have  resulted  in  two  basic 
simulation  methods. 

1.  Construction  of  a  turbulence  time  series  from  the  von 
Karman  spectra  requires  the  use  of  FFT  (Fast  Fourier 
Transform)  techniques  to  generate  the  turbulence  and 
storage  on  a  mass  storage  device  until  it  is  needed  in 
a  flight  simulation. 

2.  The  von  Karman  spectra  is  typically  approximated  to 
the  closest  integer  power  (2)  as  in  the  Dryden 
spectra,  so  that  a  set  of  difference  equations  can  be 
solved  in  real  time  during  simulation  to  generate 
turbulence  as  needed.  This  method  is  much  more 
efficient  but  provides  turbulence  that  is  less 
accurate. 

The  approximation  method  used  by  Dr.  Campbell  provides  a 
higher  order  linear  relation  which  closely  approaches  the 
-5/3  rolloff  of  the  von  Karman  spectra  as  follows. 
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The  above  approximation  was  obtained  by  truncating  at  the 
proper  point  the  continued  fraction  expansion  of  the 
binomial  function  below. 


(1+S)  =  1  +  VS  + 


(l-V)S  (1+Y)S  (2-V)S 


(N-V)S  (N+V)S 

-  +  - 
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For  a  detailed  explanation  of  this  procedure  see  "Monte 
Carlo  Turbulence  Simulation  Using  Rational  Approximations  to 
von  Karman  Spectra"  by  Dr.  C.W. Campbell (6.6). 

The  filter  function  used  in  the  transverse  case  is  derived 
in  much  the  same  manner  and  has  the  following  form. 
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The  spectra  of  these  approximations  can  be  determined  and 
compared  to  the  irrational  von  Karman  in  order  to  check  the 
validity  of  the  approximation.  Assume  that  a  stationary 
random  signal  with  a  power  spectral  density  function  4> x( w) 
is  applied  as  an  input  to  a  system  with  a  response  function 
H(w).  The  output  from  the  system  will  be  stationary  random 
signal  with  a  power  spectral  density  function  <f>y(w)  given  by 

<fry(w)  =  |  H(w)  |  2  ♦x(w) 


By  multiplying  the  filter  functions  by  their  complex 
conjugates  the  resulting  spectra  are  obtained. 


4  2 
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where  w  =  2irfal/V 
a  =  1.339 

L  =  Turbulence  scale  length 
V  =  Aircraft  airspeed 


The  following  plots  show  the  approximate  spectra  above 
plotted  versus  the  theoretical  von  Karman  spectra. 
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Figure  2.  Approximate  versus  theoretical  von  Kantian  spectra  -  longitudinal 
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Figure  3.  Approximate  versus  theoretical  von  Karman  spectra  -  transverse 


Two  computer  programs  using  this  technique  have  been 
obtained  by  Control  Dynamics  Co.  and  modified  to  run  on  a 
HP9000  computer.  The  programs  use  gausslan  white  noise  to 
drive  a  linear  filter  Implemented  as  a  set  of  difference 
equations  corresponding  to  the  von  Karman  approximation. 
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Figure  4.  Turbulence  generation  scheme  using  linear  filter  to  produce 
turbulence  having  approximate  von  Karman  spectra 


The  following  plots  show  turbulence  time  series  and  their 
calculated  PSD's  for  both  longitudinal  and  transverse  cases. 
The  turbulence  was  generated  with  gaussian  white  noise 
having  standard  deviation  of  1.0  and  mean  of  0.0.  The 
turbulence  scale  length  was  500  meters  and  vehicle  velocity 
was  100  meters/ second.  A  5.0  Hz.  sampling  frequency  was  used 
to  generate  1024  points  of  data. 
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Figure  6.  Velocity  spectra  for  longitudinal  turbulence 


7.  ATMOSPHERIC  STABILITY 


The  consideration  of  atmospheric  stability  is  important  due 
to  its  effect  on  turbulence  Intensity.  The  intensity  of 
turbulence  has  been  shown  by  Frost(6.31)  and  others  to 
Increase  with  respect  to  the  neutral  case  for  an  unstable 
atmosphere  and  to  decrease  for  a  stable  atmosphere.  The  most 
commonly  used  turbulence  spectra  (such  as  the  Dryden  and  von 
Karman)  are  valid  for  neutral  conditions. 

Atmospheric  stability  Is  measured  by  the  temperature  profile 
with  altitude.  A  decrease  in  temperature  with  altitude  is 
referred  to  as  a  lapse  rate.  An  increase  in  temperrture 
with  altitude  is  referred  to  as  an  Inversion. 

During  the  daytime  hours,  solar  radiation  heats  the  earth 
more  than  It  does  the  atmosphere.  Conduction  from  the  earth 
causes  the  air  near  the  surface  to  be  warmer  than  that 
above,  and  a  lapse  rate  results.  At  night,  with  clear  skies, 
the  earth  cools  by  radiating  heat  and  the  air  next  to  the  It 
cools  by  conduction,  thus  leading  to  1nversions(6.3) . 

Atmospheric  stability  Is  measured  by  the  tendency  of  air 
displaced  vertically  from  Its  equilibrium  condition  to 
return  to  Its  original  position.  If  the  temperature  of  a 
parcel  of  displaced  air  decreases  at  the  adiabatic  lapse 
rate  there  exists  no  restoring  force  on  the  parcel  and  the 
atmosphere  is  said  to  be  neutral.  When  a  parcel  of 
displaced  air  cools  more  quickly  than  the  adiabatic  lapse 
rate  there  is  a  positive  restoring  force  tending  to  move  the 
parcel,  which  is  now  cooler  than  the  surrounding  air,  back 
toward  its  original  position.  This  represents  a  stable 
atmosphere.  If  a  parcel  of  displaced  air  cools  more  slowly 
than  the  adiabatic  lapse  rate  there  Is  a  negative  restoring 
force  tending  to  move  the  parcel,  which  Is  now  warmer  than 
the  surrounding  air,  away  from  Its  original  position.  This 
represents  an  unstable  atmosphere. 

The  relative  occurrences  of  each  stability  category  has  not 
been  widely  studied,  although  some  data  is  available. 
Essenwanger(5.2}  compiled  a  set  of  data  taken  in  three 
locations  In  which  the  observed  atmospheric  stability  was 
recorded  for  a  one  year  period.  Numerical  tables  were 
presented  for  data  concerning  the  number  of  occurrences  of 
each  of  six  stability  categories.  The  observations  were 
grouped  by  time  of  day  (0000,0600,1200,1800  GMT  or  local), 
city  (Frankfurt,  Hahn,  Osan),  and  by  month. 

Although  this  data  Is  sufficient  to  show  general  trends 
Involved  with  the  variation  of  observed  atmospheric 
stability,  the  format  of  the  data  made  Interpretation  of  the 


data  difficult.  Control  Dynamics  has  reduced  the  data  to  bar 
chart  form  to  allow  easier  interpretation  of  trends.  The 
charts  use  the  following  conventions. 

1.  Each  chart  represents  the  occurrences  for  one  year  at 
one  time  of  day  for  one  city. 

2.  Six  different  criterion  are  established  according  to 
the  Pasquil  (5.2)  stability  criterion,  very  unstable, 
moderately  unstable,  slightly  unstable,  neutral, 
moderately  stable,  and  very  stable.  These  categories 
are  distinguished  by  texture  as  shown  on  the  charts 

1 egend. 

3.  The  number  of  occurrences  in  each  observed  stability 
category  is  represented  by  the  appropriately  shaded 
region  of  each  months  bar  chart.  The  total  height  of 
the  bar  indicates  total  observations  for  the  month, 
with  each  textured  region  indicating  the  relative 
occurrences  of  each  stability  category. 

By  studying  the  charts,  it  becomes  apparent  that  there 
exists  a  pronounced  diurnal  variation  (daily)  in  atmospheric 
stability.  In  addition  to  seasonal,  and  geographic  variations. 
The  general  patterns  are  as  follows. 

1.  In  each  city  the  atmosphere  is  neutral  to  stable  at 
midnight  for  all  seasons. 

2.  In  each  city  the  atmosphere  Is  neutral  to  unstable  at 
noon  for  all  seasons. 

3.  Each  city  shows  some  seasonal  variation  for  all  times, 
however  this  effect  is  more  pronounced  at  0600  and 
1800  when  the  effect  of  the  diurnal  variation  is  not 
as  strong. 

4.  Each  city  shows  varying  amounts  of  differences  based 
on  geographic  location. 

In  summary  the  diurnal  variation  shows  the  most  pronounced 
effect  on  stability,  with  the  seasonal  and  geographic 
effects  of  lesser  Importance  respectively. 

The  following  figure  obtained  from  Frost(6.31)  shows 
measured  data  for  stable,  unstable,  and  neutral 
atmosphere's.  This  plot  Indicates  the  changes  In  the  shape 
of  the  observed  spec trim. 
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9.  Measured  data  for  stable,  neutral ,  and  unstable  atmospheres  vs  the 
theoretical  von  Karman  and  Dryden  spectra  (Wave  number  K  *  w/v 
normalized  by  Jhe  peak  scale  Km  of  the  vertical  spectra.) 

Note:  f$(f)/o  where  a  *  standard  deviation  f  *  frequency. 

As  can  be  seen,  the  spectra  of  the  stable,  neutral,  and 
unstable  data  are  the  same  and  correspond  to  the  von  Karman 
spectra  at  frequencies  above  the  -5/3  slope  break  point.  The 
value  of  this  break  point  is  dependent  upon  turbulence  scale 
lengths  and  vehicle  or  reference  velocity.  At  lower 
frequencies  the  data  diverges  somewhat  from  the  neutral 
case.  The  figure  clearly  shows  the  effect  of  a  stable 
atmosphere  decreasing  turbulence  intensity  and  an  unstable 
atmosphere  Increasing  turbulence  intensity.  The  neutral 
atmosphere  corresponds  to  the  theoretical  von  Karman  and 
Dryden  spectra. 


In  order  v  simulate  the  effects  of  stable  and  unstable 
atmospheres  while  using  models  with  spectra  corresponding  to 
a  neutral  atmosphere,  some  modifications  need  to  be  made. 

One  such  modification  that  is  particularly  applicable  to 
Monte  Carlo  simulations  is  the  addition  of  a  low  pass  filter 
path  so  that  the  desired  lower  frequencies  could  be  scaled 
by  a  chosen  gain  factor  as  below 
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Figure  10.  Method  for  scaling  turbulence  spectra  for  effects  of  atmospheric 
stabil ity. 


The  design  of  this  low  pass  filter  would  have  to  be 
carefully  chosen  so  as  not  to  overly  distort  the  high 
frequency  part  of  the  turbulence  spectrum. 


Another  method  of  accounting  for  the  effects  of  atmospheric 
stability  on  turbulence  intensity  involves  the  use  of  the 
Kaimal  turbulence  spectra  as  presented  in  Frost(6.31).  The 
Kaimal  spectra  consists  of  an  empirical  formula  as  follows. 


*(f)  = 


i  [0.164  (n/no)  ] 

5/3 

f  [1  +  0.164(n/no)  ] 


f  =  cyclic  frequency 
a  =  standard  deviation 
n  =  the  reduced  frequency  =  fz/V 
z  =  height 

V  *  reference  velocity 


The  values  of  "no" 
conditions  are 


recommended  by  Frost  for  neutral 


noj  *  0.0144 
n02  =  0.0265 
n03  -  0.0962. 

By  changing  these  values  the  shape  of  the  Kaimal  spectra  may 
be  taylored  to  different  atmospheric  stabilities. 

Other  authors  have  suggested  a  much  simpler  although  less 
accurate  method  of  accounting  for  stability  effects  on 
turbulence. (5.2)  The  method  involves  the  scaling  of  the 
wind's  standard  deviation  as  a  function  of  stability.  Since 
the  turbulence  spectra  is  directly  proportional  to  the 
standard  deviation,  this  modification  causes  the  spectra  to 
be  shifted  up  or  down  with  respect  to  the  neutral  case.  This 
results  in  a  spectra  which  becomes  more  accurate  at  low 
frequencies  but  less  accurate  at  high  frequencies.  Whether 
this  method  is  acceptable  to  a  given  problem  can  only  be 
determined  through  a  knowledge  of  the  frequency  ranges  that 
are  of  importance  to  the  problem  being  considered. 
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9.  SPATIAL  MODELING 

The  topic  of  spatial  modeling  of  wind  is  one  in  which 
significantly  less  work  has  been  performed.  The  primary 
topics  of  concern  for  spatial  turbulence  models  are  the 
coherence  of  turbulence  along  with  Taylor's  hypothesis. 

Over  somewhat  short  distances  (<  100-200  m)  turbulence  can 
be  shown  to  be  coherent  to  some  degree  in  each  axis.  This 
coherence  is  typically  represented  as  an  exponentially 
decaying  function  such  that  the  turbulence  becomes  more 
rapidly  less  coherent  as  the  separation  distance  increases. 
Typical  coherence  curves  as  taken  from  Armendariz(3.1)  are 
shown  in  the  following  figures. 


4  \  * 
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Fiqure  11.  Zero-Lag  horizontal  wind  component  correlation  at  a  height  of  1.5 
meter  under  unstable  conditions  as  function  of  distance  separa¬ 
tion  for  mean  wind  parallel  and  normal  to  sensor  line. 
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Figure  13.  Zero-Lag  horizontal  wind  component  correlation  at  a  height  of  16 
meter  under  unstable  conditions  as  function  of  distance  separa¬ 
tion  for  mean  wind  parallel  and  normal  to  sensor  line. 

The  rate  at  which  the  coherence  decays  Is  a  function  of 
separation  distance,  frequency,  mean  wind  speed,  and  some 
decay  constant  as  shown  below. 


-{awx/v) 

coh  =  e 

where  a  *  decay  constant 

w  *  turbulence  frequency 
x  *  separation  distance 
v  *  reference  velocity 

Each  of  these  terms  Is  relatively  straightforward  except  for 
the  decay  constant  "a".  While  all  authors  include  this  term, 
few  If  any  agree  on  its  value.  It  seems  to  be  somewhat 
arbitrarily  chosen  based  on  experience  and  data  available. 


V7 

s'  1 
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This  expression  for  coherence  allows  the  following  generalizations 

1.  Strong  winds  show  greater  coherence  than  light  winds. 

2.  Low  frequency  eddies  show  greater  coherence  than  high 
frequencies. 

3.  Coherence  decreases  Increasingly  fast  with  separation 
distance. 

4.  Coherence  is  greater  In  an  unstable  atmosphere. 

The  modeling  of  coherence  generally  Incorporates  the  use  of 
many  Independent  noise  sources  which  are  filtered  and  added 
before  driving  some  turbulence  filter  function.  The 
appropriate  phase  lagging  of  these  input  filters  based  on 
separation  distance,  reference  velocity,  and  frequency 
produces  correlated  turbulent  time  series  for  a  given 
separation  distance.  The  form  of  the  inter-level  coherence 
model  developed  by  Frost(6.31)  is  shown  below. 


Figure  14.  Inter-level  coherence  model 

While  this  type  of  turbulence  will  have  a  significant  Impact 
on  the  study  of  structures  or  on  the  dynamics  of  large 
vehicles.  It  would  seem  that  of  greater  Importance  to 
rocketry  problems  Is  a  discussion  of  Taylor's  hypothesis. 

Taylor's  hypothesis  Involves  the  assumption  that  turbulence 
Is  frozen  In  the  mean  wind,  because  the  turbulence  decays 
slowly  with  respect  to  the  speed  at  which  the  observer 
passes  through  It.  In  other  words,  we  might  obtain  about  the 
same  data  by  recording  In  time  at  a  fixed  location  as  we 
would  by  rapidly  moving  upstream  against  the  mean  wind. 


To  Illustrate  these  Ideas  consider  the  record  of  an 
atmospheric  variable  A,  observed  at  a  fixed  point  as 
depicted  In  the  following  figure.  On  the  assumption  that 
turbulent  eddies  are  carried  by  the  mean  wind  "U"  and  change 
only  slowly  as  they  move,  we  could  assume  that  the  pattern 
In  the  figure  is  produced  by  having  a  spatial  pattern  in  the 
"x"  direction  moving  past  the  observer  with  speed  "U“ 
without  change  of  shape.  In  other  words,  the  spatial  pattern 
could  be  depicted  exactly  as  in  the  figure  with  the  axis 
relabeled  as  "x=Ut".  This  notion  is  called  the  frozen  wave 
hypothesis  or  Taylor's  hypothesis(7.5). 


t 


Figure  15.  Illustration  of  Taylor’s  hypothesis 


As  is  evident  from  the  definition  above,  in  order  that 
Taylor's  hypothesis  be  valid,  it  is  necessary  that  the 
turbulence  be  stationary  in  time  and  homogeneous  in  the  x 
direction.  These  conditions  are  often  satisfied  near  the 
boundary  layer  when  the  terrain  is  consistent  in  roughness. 

Work  on  the  experimental  validation  of  Taylor's  hypothesis 
has  been  performed  by  Frost(5.9)  at  the  Marshall  Space 
Flight  Center  boundary  layer  tower  array.  In  this 
experiment  turbulence  time  series  were  taken  at  several 
towers  aligned  with  the  wind.  This  actual  data  was  then 
compared  to  data  which  was  taken  at  the  upwind  tower  and 
scaled  by  the  mean  wind  speed  to  produce  a  spatial 
turbulence  series.  According  to  Frost(5.9),  this  frozen  and 


Lag  Time,  t  (sec) 


H 

Figure  16.  Illustration  of  validity  of  Taylor's  hypothesis 
(•  spatial  correlation  converted  to_t1me 
correlation  with  the  relationship  z=wt). 


scaled  turbulence  gave  excellent  agreement  to  what  was 
actually  measured  at  various  points  downwind.  The  following 
figure  from  Frost  illustrates  the  agreement  between  actual 
data  and  data  predicted  by  Taylor's  hypothesis. 

Thus  It  would  seem  that  for  terrain  that  Is  homogeneous  in 
nature,  Taylor's  hypothesis  provides  a  useful  tool  for 
obtaining  a  spatial  wind  history  from  a  wind  measurement. 

A  second  point  to  be  considered  with  Taylor's  hypothesis 
Involves  the  simulation  of  a  turbulent  time  series.  For  the 
case  In  which  a  vehicle  Is  moving  through  a  turbulent  field 
rapidly,  the  turbulence  can  be  considered  frozen.  Then 


scaling  the  time  series  by  the  vehicle  airspeed  produces  a 
spatial  series.  Therefore  for  this  type  of  situation  a 
turbulence  time  series  can  be  considered  as  spatial 
turbulence  simply  by  scaling. 

In  summary,  the  consideration  of  correlation  and  Taylor's 
hypothesis  are  Important  for  cases  such  as  structure 
loading,  where  the  separation  distance  between  points  of 
Interest  are  relatively  small.  To  expect  a  high  correlation 
of  turbulence  2  to  3  kilometers  from  a  measurement  site  is 
however  unrealistic.  According  to  Armendariz(3.1),  the 
coherence  decreases  to  almost  negligible  values  with 
horizontal  separation  greater  than  200  meters  even  under 
unstable  conditions,  except  for  the  long  wavelength  portion 
(>3000m).  This  suggests  that  one  can  only  extrapolate  the 
long  wave  portion  of  the  wind  energy  beyond  a  few  seconds 
and  for  greater  distance  than  perhaps  200  meters.  The  short 
wave  or  turbulent  portion  appears  to  be  quite  intractable. 

It  becomes  apparent  that  to  predict  significantly  the  wind 
conditions  In  the  target  area  requires  the  use  of  remote 
measurements.  Perhaps  the  most  promising  method  of 
obtaining  remote  measurements  without  the  use  of  downrange 
Instruments  is  the  utilization  of  LIDAR.  Frost(5.8)  has 
reported  on  the  use  of  HOAR  In  conjunction  with  an 
Instrumented  aircraft  In  the  NASA  gust  gradient  program. 

The  LIDAR  was  used  to  take  remote  measurements  of 
turbulence  along  the  ILS  approach  path  and  was  later 
compared  to  the  data  taken  aboard  the  aircraft  flying  the 
approach.  Good  agreement  was  achieved  between  the  LIDAR  and 
the  aircraft  measurements. 


10.  FUTURE  WORK 


The  next  step  would  probably  be  the  development  of  accurate  mean 
wind  and  turbulence  models,  suitable  for  implementation  into  a 
6 DOF  guided  missile  simulation.  Due  to  the  variety  of  simulations, 
the  model  would  be  Incorporated  into  a  designated  system.  The  model 
would  be  capable  of  handling  stable,  neutral,  or  unstable  atmospheric 
conditions,  Gaussian  or  non-Gausslan  distributed  turbulence,  general 
terrain  roughness  characterization,  gust  gradient  effects,  and 
coherent  turbulence. 

Models  can  be  developed  for  tree  line  created  turbulence  or  other 
major  obstacles  and  verified  with  a  three  dimensional  water  tunnel 
flow  simulation  capability  such  as  can  be  performed  by  FWG. 

To  access  accurately  the  effects  of  variable  wind,  the  standard 
equations  of  motion  used  to  describe  missile  dynamics  need  to  be 
expanded  since  they  typically  assume  either  no  wind  or  a  constant 
mean  wind.  The  additional  terms  would  include  both  temporal  and 
spatial  gradients  of  wind  (Bowles  and  Frost  7.3).  The  relative 
Importance  of  these  terms  and  related  changes  to  input  parameters 
for  the  aerodynamic  coefficients  need  to  be  assessed  on  a  missile 
by  missile  basis. 

The  Incorporation  of  accurate  wind  turbulence  models  allows  for  a 
more  realistic  assessment  of  dispersion  of  dumb  rounds.  One 
Interesting  possibility  would  be  to  try  to  optimize  the  design  of 
the  round  to  minimize  the  dispersion. 

Since  wind  prediction  based  on  correlation  has  been  shown  to  be 
unrealistic  beyond  one  or  two  hundred  meters,  the  alternative  of 
using  LIDAR  to  measure  data  along  the  anticipated  trajectory 
should  be  examined.  Simulations  would  be  used  in  the  first  phase 
to  assess  potential  gains,  followed  by  a  subsequent  phase  that 
would  Involve  actual  use  of  LIDAR  with  the  RAKE  or  some  similar 
system  that  would  allow  for  measurements  to  be  taken  at  other  than 
the  launch  site. 


11.  SUMMARY 


Control  Dynamics  has  performed  as  extensive  survey  of  the 
literature  available  and  pertinent  to  the  problems  associated 
with  wind  effects  on  free  flight  rockets.  In  addition  to 
surveying  the  available  literature,  meetings  and  phone  contacts 
have  been  held  with  Individuals  who  are  actively  involved  with 
research  in  the  field  of  atmospheric  science.  The  primary 
Intent  of  these  surveys  and  contacts  has  been  to  determine  what 
type  of  work  Is  being  performed  that  has  an  Impact  on  design 
and  simulation  of  small  missiles.  It  was  also  desired  to 
establish  areas  which  seem  to  be  particularly  Important  and 
which  warrant  further  study. 

Several  areas  which  seem  to  be  of  Interest  Include  the  use 
of  the  Space  Shuttle  Turbulence  Tapes  (SSTT),  the  approximation 
of  the  von  Karman  spectra  with  a  linear  expression,  the  effects 
of  atmospheric  stability  on  turbulence  intensity,  mean  wind 
data,  and  several  Issues  concerning  spatial  modeling. 

The  Space  Shuttle  Turbulence  Tapes  provide  a  quick  and  rela¬ 
tively  simple  means  of  Incorporating  high  quality  turbulence 
data  into  current  simulations.  Reservations  concerning  the 
lack  of  high  frequency  content,  mass  storage  problems,  and 
necessity  of  always  using  the  same  turbulence  would  seem  to 
Indicate  that  the  SSTT  are  primarily  valid  for  large  and 
slower  vehicles. 

The  von  Karman  approximation  method  provides  a  compact,  effi¬ 
cient,  means  of  generating  turbulence.  The  code  Is  efficient 
enough  to  generate  turbulence  on  an  as  needed  basis  Instead 
of  storing  large  blocks  of  data  in  mass  storage.  This  pro¬ 
vides  an  added  benefit  of  generating  Independent  sets  of 
turbulence  by  changing  the  seeds  on  the  random  number  genera¬ 
tor.  As  shown  In  the  plots  in  the  body  of  this  report,  the 
turbulence  provides  a  very  close  fit  to  the  von  Karman  at  all 
but  very  high  frequencies. 

A  discussion  of  atmospheric  stability  along  with  data  showing 
the  type  of  changes  that  occurs  on  a  daily,  seasonal,  and 
geographical  basis  has  been  Included.  This  data  has  been 
presented  in  a  bar  chart  form  Instead  of  the  original  tabular 
form  In  order  to  provide  Increased  clarity  and  ease  of  Inter¬ 
pretation.  The  data  shows  that  the  changes  in  atmospheric 
stability  are  sufficient  to  warrant  Including  In  a  realistic 
wind  turbulent  model.  The  effects  of  stability  are  typically 
Ignored  In  most  turbulence  models,  but  several  methods  exist 
which  would  allow  the  accurate  modeling  of  stability  effects 
on  turbulence  intensity.  These  Include  scaling  the  spectra 


by  a  gain,  low  pass  filtering  the  turbulence,  or  the  use 
of  an  alternate  spectra  such  as  the  Kaimal  which  includes 
stability  terms. 


Data  for  mean  winds  and  standard  deviations  at  35  sites  in 
the  northern  hemisphere  Is  available  in  tabular  form.  This 
data  was  presented  in  bar  chart  form  for  3  sites  in  order  to 
show  a  representative  sample  of  the  available  data. 


Various  topics  have  been  discussed  regarding  spatial  model¬ 
ing  of  wind.  These  Include  the  Inclusion  of  coherence  and 
Taylor's  hypothesis  Into  a  model  and  the  usefulness  of 
these  techniques.  The  use  of  coherence  modeling  and  Taylor's 
have  been  shown  to  be  valid  and  important  for  situations 
concerned  with  the  modeling  of  wind  over  relatively  short 
d i stances (<200m) .  It  becomes  evident  however,  that  at  large 
distances  such  as  occur  in  rocketry  problems  (2-3km)  these 
techniques  are  no  longer  valid.  The  best  that  can  be 
expected  in  wind  prediction  over  long  ranges  is  to 
characterize  the  general  statistics  of  the  wind. 


The  local  measurements  allow  the  rough  prediction  of  mean 
wind  and  turbulence  intensity  along  with  the  basic  spectra 
that  is  being  observed.  This  of  course  assumes  homogeneous 
terrain.  If  the  terrain  has  large  irregularities  between 
launch  and  target  sites,  then  it  may  be  Impossible  to 
predict  accurately  the  effects  of  wind  on  the  rocket 
downrange.  The  basic  problems  Involving  long  range  wind 
predictions  are  the  following. 

1.  Turbulence  degrades  excessively  as  it  moves  downrange 
due  to  the  long  time  it  takes  to  be  blown  downwind. 


2. 


Terrain  is  more  likely  to  be  Inhomogeneous  over  long 
distances. 


3.  Wind  is  unlikely  to  be  aligned  along  the  path  of  the 
missile.  In  this  case  the  measurement  may  be  biased  by 
terrain  near  the  launch  site.  This  terrain  Is  likely 
to  be  different  from  the  terrain  along  the  flight  path. 


Additional  Improvements  In  wind  measurements  could  be  made 
by  using  advanced  instrumentation  such  as  HOAR  to  give 
directly  measured  downrange  wind  measurements.  Additional 
sophistication  In  modeling  could  be  obtained  through  the  use 
of  advanced  numerical  models  In  combination  with  water  tunnel 
simulation  of  complex  terrain. 
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13.  APPENDIX 

This  appendix  contains  listings  to  the  computer  programs 
used  to  produce  the  turbulence  time  series  by  approximation 
of  the  von  Karman  spectra.  These  programs  were  written  and 
supplied  by  Dr.  Warren  Campbell  at  Marshall  Space  Flight 
Center  and  were  modified  by  Control  Dynamics  to  run  on  a 
HP9040. 
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****************  LONGITUDINAL  ***************** 

THIS  PROGRAM  GENERATES  “TURBULENCE"  USING  A  2-3  RATIONAL 
APPROXIMATION  TO  THE  VON  KARMAN  SPECTRUM. 

PROGRAMMER:  OR.  WARREN  CAMPBELL  <  NASA  MARSHALL  SPACE  FLIGHT  CENTER) 

REVISOR:  DAVID  EDGEMON  <  CONTROL  DYNAMICS  CO.  ) 

C - 

c  DIMENSION  IPARC5) 

REAL*8  YO ,  Y1 ,  Y2 , YN, XG ,  X 1 , X2, X3, ELSO, ELS1 , ELS2, ELS3, RO, Rl , R2, R3 
8.  ,DT,PI,EX1  ,EX2,EX3,EX12,EX13,EX23,EX123 
REAL+8  PI ,P2,P3,Z1 ,22, A, 8,0,0, XK1 , SUM, SUMSQ, YBAR, Y2BAR 
INTEGER+4  ITER, NPTS 
integer  irGw,ix,iy 
real*8  sig, xb, xdp, xp 

real +8  avk ,alpha,el,sl , sigma, v, time, yp 
c 

integer  nrec,icol,j 
real  outv< 2  ) 
c 

COMMON/CFF/  SIGMA, V, ITER , AVK , EL , DT , PI 
COMMON/EXPBL/  EX1 , EX2, EX3, EXt 2, EXt 3, EX23, EX t 23 
c 

open<  unit=l , form- 'unformatted ' , f i le=  'plotl . dat ' , access5 'direct 
&, recl«4  ) 

open<  unit =2, f i le=  'plot 2 . dat '  ) 
open<  unit =3, f i le=  'time . dat ' > 
c 

icol=2 

irow=0 

c 

write<6,*)  'enter  standard  deviation  for  random  numbers' 
read<5,*)  sig 

write<6,*)  'enter  mean  of  random  numbers' 
read<5,*)  xb 

write<6,*)  'enter  seed  for  random  numbers' 
read<5,*)  ix 
c 

c  CALL  RMPAR(IPAR) 

C - 

c  LUPRT  =»  I PAR<  1  ) 

C - 

C 

C  SET  UP  CONSTANTS 
C 

AVK  =  1.339 
UIRITE<  6 , 5999  > 

5999  F0RMAT<  '  ENTER  HUMBER  OF  POINTS') 

READ<5,*)  NPTS 
WRITE<  6, 7998  ) 

7993  F0RM*T< '  ENTER  DT  IN  DOUBLE  PRECISION') 

READC  5, *  )  DT 
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PI  =  4 . +DATANC 1 .DO) 

ITER  -  0 
SUM  =  O.DO 
SUMSQ  =  0.0D0 
WRI TEC  6, 7999  ) 

7999  FORMATC  '  ENTER  DESIRED  SIGMA') 

READC5,*)  SIGMA 
SI  =  SIGMA 

D  WRITEC6,8999)  PI 

8999  FORMATC  '  PI  =  ' , D1 5 . 1 0  ) 

C 

C  INITIALIZE  THE  X'S  AND  Y'S 
C 

CALL  INITalCYQ,Y1,Y2,X0,X1,X2,X3,TIME,  ITER,  sig,xb,xdp,xp,  ix,  iy) 
D  WRITEC  6, 8998  )  TIME 

8998  FORMATC  '  TIME  = ' ,  FI  0 , 5  ) 

C 

C  GET  SIGMA,  V,  AND  TURBULENT  LENGTH  SCALE 
C 

100  CONTINUE 

CALL  SVLCDT, SIGMA, V, EL) 

ALPHA  *  V*DT/EL 

SIGMA  *  S17C ,8762*ALPHA**0. 49626) 

D  WRITEC 6 , 8997  ) 

8997  FORMAT<  '  SVL  '  ) 

C 

C  CALCULATE  POLES  AND  ZEROES 
C 

CALL  PANDZ<P1 ,P2,P3,Z1,Z2) 

D  WRI TEC 6, 8996 >  PI , P2, P3, 21 , Z2 

8996  FORMATC  '  POLES= ' , 3< 2X, D1 5 . 9 >/  '  ZEROES* ' , 3<  2X, D1 5 , 9  > ) 

C 

C  CALCULATE  COEFFICIENTS 
C 

CALL  COEFFCP1 ,P2,P3,Z1 , Z2 , A, B , C, D , XK 1  ) 

D  WRITEC6,8995)  A,B,C,D,XK1 

8995  FORMATC'  COEFF  A, B, C, D, XK1  = ' , 5< 2X , El  5 . 9 >  ) 

C 

C  CALCULATE  RHS  DIFFERENCE  EQUATION  COEFFICIENTS 
C 

CALL  RIC PI , P2,P3,A,B,C,D,R0,R1 , R2, R3 ) 

D  WR I  TEC  6 , 8994  )  R0,R1,R2,R3 

8994  FORMATC'  RI= ' , 4C 2X, D1 5 . 9 > > 

C 

C  CALCULATE  DIFFERENCE  EQN  LEFT  SIDE  COEFFICIENTS 
C 

CALL  LICP1 , P2, P3, A, B, C, D, ELS0,ELS1 , ELS2, ELS3 ) 

D  WR I  TEC  6, 8993  )  ELSO, ELS  1 , ELS2 , ELS3 

8993  FORMATC'  LI* ' , 4C  2X, D1 5 . 9  >  > 

C 

C  ITERATE 
C 

YN  *  ELS1 *Y2  -  ELS2*Y1  +  ELS3*Y0  +  R0*X3  -  R1*X2  +  R2*X1 
«.  -  R3*X0 

ITER  *  ITER  +  1 
SUM  *  SUM  +  YN 
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SUMSQ  =  SUMSQ  +  YN+YN 
YBAR  =  SUM/ITER 

Y2BAR  =  OSQRT<  SUMSQ/ITER  -YBAR**2 > 

D  URITE<  6, 8900  >  Y2,  Y1  ,  Y0,  X3,  X2,  X 1  ,  XO 

8900  FORMAT<  '  Y2 , Y 1 , Y  0= ' , 3<  2X , D 1 5 . 9  V '  X3 , X2 , X 1 , X 0=  ' , 4< 2X , D 1 5 . 9 >  ) 
YO  «  Y 1 
Y1  *  Y2 
Y2  =  YH 
XO  -  XI 
XI  «  X2 
X2  =>  X3 

call  random  < sig, xb , xdp , xp , ix , i y  ) 

X3  *  xdp 

D  WRI TE<  6 , 9998  )  X3 

9998  FORMAT<  '  RANDOM  NUMBER®  ' ,  FI  0 . 6  > 

TIME  =  TIME  +  DT 
YP  =  YN 

IJR1  TE<  6 , 9999  )  TIME,  YP 
c 

c  write  data  to  plot  files 

c 

outv( 1  )=time 
outv<  2  }®yp 
do  1 234  j® 1 , icol 
nrec=nrec+ I 

write( 1 ,rec»nrec)  outv< j > 

1234  continue 

c 

c  write  data  to  file  for  psd  program 

c 

write<3, '<e13.7>'>  yp 
iron® irow+1 


9999  FORMAT<  '  TIME®  ' , El  2 . 6,  '  Y=',E15.9> 

IFCITER  .LT.  NPTS )  GO  TO  100 
WRITE<  6, 9980 )  NPTS, YBAR, Y2BAR 
9980  F0RMAT<  '  NPTS®', 18,'  YBAR= ' ,  D1  5 . 9,  '  Y2BAR®  ' , D1 5 . 9  > 
c 

write<2, '< 2i5  )  '  )  icol,irow 
write<2,  '< aS >  '  )  'time' 
write<2,  '<a8>'  )  'turb' 
c 

c  IFCLUPRT  .EQ.  8)  ENDFILE  8 

c 

close<unit=1 ) 
close<  unit  =  2 > 
c 

STOP 

END 

c 

c 

SUBROUTINE  INITal< Y0, Y1 , Y2 , X0, X 1 , X2, X3, T IME, ITER 

t, sig, xb, xdp, xp , ix , iy  > 

REAL+8  Y0,Y1 ,Y2,X0,X1 ,X2,X3,time 

INTEGER*4  ITER 
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c 

c 


c 

c 


c 

c 


c 


real*8  s ig, xb , xdp , xp 
integer  ix,  iy 
XO  =  0  .  DO 
XI  =  0 . DO 
X2  =  0 . DO 

call  random  <  s  i  g ,  xb ,  xdp ,  xp ..  i  x ,  i  y  > 
X3  *  xdp 
Y0  *  0 . DO 
Y 1  =  0 . DO 
Y2  =  0 . DO 
TIME  =  0. 

ITER  =  0 

RETURN 

END 


SUBROUTINE  SVL<  DT, SIGMA, V, EL  ) 
real *8  dt,  sigma.,  el ,  si  ,  v 
SI  =  SIGMA 
V  =  100. 

EL  =  500. 

RETURN 

END 


SUBROUTINE  PANDZ< P 1 , P2, P3, Zl , Z2 > 

COMMON/CFF/  SIGMA, V, ITER, AVK, EL, DT, PI 
REAL*8  PI , P2, P3, 21 , Z2, C, PI , DT 

real *8  avk , el , pi p, p2p, p3p, sigma, v,  zl p, z2p 
INTEGER*4  ITER 

DATA  Z1P/-1 .46821 076/, Z2P/-5 . 388932 1 43/ , P 1 P/- 1 . 02025059/ 
*, P2P/-1 . 67661571/, P3P/-8. 1 0313370/ 

C  «  V/<  AVK*EL  ) 

PI  =  PI P*C 
P2  =  P2P*C 
P3  =  P3P+C 
Zl  =  Z1P*C 
Z2  =  Z2P*C 
RETURN 
END 


SUBROUTINE  COEFFC PI , P2, P3, Zl , Z2, A, B, C , D, XK1  ) 
COMMON/CFF/  SIGMA, V, I  TER , AVK, EL , DT , P I 
REAL*8  PI ,P2,P3,Z1 ,Z2,A,B,C,D,XK1 ,DT,PI 
real *8  avk , el , sigma, v 
INTEGERM  ITER 

XK1  -  SIGMA*DSQRT<  EL/<  V*PI  >  ) 

XK1  *  XK1*<-P1*P2*P3>/<Z1*Z2> 

A  *  XK1  *<  PI  -Zl  >*<  PI  -Z2  >/<  PI *<  PI  -P2  )*<  PI  -P3  >  > 
B  «  XK1 *<  P2-Z 1 >*<  P2-Z2  >/<  P2  *<  P2-P 1 >*( P2-P3  )  > 
C  *  XK1*<P3-Z1  )*< P3-Z2 >/< P3*< P3-P1  >*<P3-P2>) 
D  -  -XK1*Z1*Z2/<P1*P2*P3  ) 

RETURN 

END 
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c 

SUBROUTINE  RIC  PI , P2,P3,A,B,C,D,R0,R1 , R2, R3  > 

REAL*8  DT, PI, PI  ,P2,  P3,  A,  B,  C,  D,  fi 0 ,  R 1  ,R2,R3,EX1  ,  EX2 ,  EX3 ,  EX  1  2 ..  EX  1  3 
&  , EX23, EX1 23 
INTEGER*4  ITER 

real*8  avk , el , sigma,  v 
COMMON/CFF/  SIGMA, V, ITER, AVK, EL, DT, PI 
COMMON/EXPBL/  EX1 , EX2, EX3, EX  1 2, EX  1 3, EX23, EXJ 23 
EX1  =  DEXP<P1*DT> 

D  UR  I TE<  6 , 9999  )  EX1 
9999  FORMATC '  EX1  =  ',D15.9> 

EX2  =  DEXPCP2*0T> 

D  UR  I TEC  6 , 9998  )  EX2 

9998  FORMATC  '  EX2=',D15.9) 

EX3  =  DEXP<  P3*DT  ) 

D  URITE<6,999?>  EX3 

999?  FORMATC '  EX3=',D15,9> 

EX1 2=  DEXPC  <P1+P2>*QT) 

D  URITEC6,9996  )  EX12 

9996  FORMATC'  EX1 2= ' , D1 5 . 9 > 

EX1  3=  DEXPC  C  P 1 +P3  >*DT  > 

D  URITEC.6,9995  )  EX13 

9995  FORMATC'  EX1 3=  ' , D 1 5 . 9  ) 

EX23=  DEXPC  CP2  +  P3>*DT> 

D  UR  I  TECS, 9994  )  EX23 

9994  FORMATC'  EX23  = ' , D 1 5 . 9 > 

EX1 23=  DEXPCCP1+P2+P3)*DT> 

D  WR I  TEC  6, 9993)  EX 123 

9993  FORMATC'  EX1 23= ' , D1 5 . 9 > 

RQ  =  A  +  B  +  C  +  D 
D  URITEC  6, 9989  )  R0 

9989  FORMATC'  R0=',D15.9> 

R1  =  A*C 1 . +EX2+EX3 )+B*< 1 . +EX 1 +EX3  >+C*C 1 . +EX 1 +EX2  >+ 

&  D*C  EX1 +EX2+EX3  ) 

R2  *  A*C  EX2+EX3+EX23  )+B*<  EX1 +EX3+EX! 3)+C*CEX1 +EX2+EX1 2  >+ 

*.  D*C  EX  1  2+EX 1  3+ EX 23  ) 

R3  =  A*EX23+B*EX1 3+C*EX1 2+D*EX1 23 
D  UR  I  TEC  6 , 9992  >  R0,R1,R2,R3 

9992  FORMATC'  R0=',D15.9,'  R1=',D15.9,'  R2=',D15.9,'  R3=',D15.9) 

RETURN 
END  . 
c 
c 

SUBROUTINE  L IC PI , P2 , P3 , A, B, C , D , ELS 0, ELS  1 , ELS2, ELS3  ) 

REAL*8  EX1 , EX2, EX3, EX1 2, EX1 3, EX23 , EX 1 23 , P I , DT 
REAL*8  PI , P2, P3, A, B, C, D, ELSO, ELS1 , ELS 2, ELS 3 
real*8  avk , el , sigma, v 
COMMON/CFF/  SIGMA, V, I  TER , AVK, EL , DT, PI 
COMMON/EXPBL/  EX1 , EX2, EX3 , EX1 2, EX  1 3, EX23, EX  1 23 
INTEGER*4  ITER 
ELSO  -  1 .DO 

ELS1  *  EX1  +  EX 2  +  EX 3 
ELS2  -  EX1 2  +  EX  1 3  +  EX23 
ELS3  -  EX  1 23 
RETURN 
END 
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c 

c 

BLOCK  DATA  EX 

REAL+8  DT,PI,EX1 , EX2, EX3, EX1 2, EX1 3, EX23, EX1 23 
INTEGER*4  ITER 

real*8  avk , e 1 , s igma, v 
COMMON/CFFX  SIGMA, V, ITER, AVK, EL, DT, PI 
COMMON/EXPBL/  EX  1 , EX2, EX3, EX  t  2 , EX  1 3 , EX23, EX  1 23 
END 
c 
c 
c 

subroutine  random  <sig, xb, xdp, xp, ix, iy > 
real*8  xp, a, xdp, sig, xb 
integer  iy,ix,j 
c 

a=  0 . 0 

do  30  j=1 ,12 
c 

iy=ix*65539 
if<iy>  10,20,20 
10  iy=iy+21474S3647+1 

20  xp=iy 

xp=xp/21 47483647 . dO 
ix=iy 

30  a=a+xp 

xdp=< a-6 . 0  )*s ig  +  xb 
c 

return 


o  o  cn  n  n  n  nonon  n  n  n  n  ooooonooo 
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PROGRAM  DTRNS 


****************  TRANSVERSE  ***************** 

THIS  PROGRAM  GENERATES  ■'TURBULENCE'*  USING  A  3-4  RATIONAL 
APPROXIMATION  TO  THE  VON  KARMAN  TRANSVERSE  SPECTRUM. 

PROGRAMMER:  DR.  WARREN  CAMPBELL  CNASA  MARSHALL  SPACE  FLIGHT  CENTER) 

REVISOR:  DAVID  EDGEMON  <  CONTROL  DYNAMICS  CO.) 


DIMENSION  I PAR (5) 

REAL*8  YO ,  Y 1  , Y2 , Y3 , YN, XO , XI  , X2 , X3 , X4 , ELS  0 , ELS1  ,  ELS2 ,  ELS3 ,  ELS4 
8,  ,R0,R1 ,R2,R3,R4 

&  , DT , P I , EX1  ,  EX2,  EX3 ,  EX4 , EX 12, EX 13, EX  14,  EX23 ,  EX24 ,  EX34  ,  EX1  23 ,  EX  1  24 
&  ,EX134,EX234,X1234 

REAL*8  PI ,P2,P3,P4,Z1 ,  Z2,  Z3, A , B, C, D , F, XK 1 , SUM , SUMSQ , YBAR ,  Y2BAR 
real  *8  avk  ,  alpha,  el, si  ..  sigma,  v,time;yp 
I NTEGER*4  I  TER , NPTS 

real*8  sig, xb, xdp,xp 
integer  ix, iy, irow, icol 

integer  nrec,j 
real  outv<2) 

COMMOH/CFF/  SIGMA, V, ITER , AVK, EL , DT , P I 

COMMON/EXPBL/  EX1 , EX2 , EX3 , EX4 , EX1 2, EX1 3 , EX1 4 , EX23 , EX24 
&, EX34, EX123,EXt  24 , EX1 34, EX234 , X 1 234 
CALL  RMPARC IPAR ) 


LUPRT  =  IPARC 1  ) 


open<  uni t  =  1 , f orm=  'un formatted', file='plot1.dat',access='direct' 
&, recl=4 ) 

open<  unit=2,file='plot2.dat'  > 
open<unit=3,file='time.dat'  ) 

i row*  0 
icol=»2 

write<6,*>  'enter  standard  deviation  for  random  numbers' 
read<5,*)  sig 

write<6,*>  'enter  mean  of  random  numbers' 
read<5,*>  xb 

write<6,*>  'enter  seed  for  random  numbers' 
read<5,*>  ix 

WR I TE<  6 , 5999  ) 

999  FORMAT< '  ENTER  NPOINTS') 

READ<  5, *  )  NPTS 


SET  UP  CONSTANTS 
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AVK  =  1 .339 
WRITEC  6.,  7998  ) 

7998  FORMATC  '  ENTER  DT  IN  DOUBLE  PRECISION') 

READC5,*)  DT 

PI  =  4.*DATANC 1 .DO) 

ITER  =  0 
SUM  *  0  .DO 
SUMSQ  =  O.ODO 
WRITEC6,7999  ) 

7999  FORMATC  '  ENTER  SIGMA') 

READC  5 ,  *  >  SIGMA 

SI  =  SIGMA 

D  WR I  TEC  6  ..8999  )  PI 

8999  FORMAT C  '  PI*',Q15.10> 

C 

C  INITIALIZE  THE  X'S  AND  Y'S 
C 

CALL  IHITalCYO,Y1  ,Y2,  Y3..XO,.  XI  , X2,X3,X4,  TIME,  ITER) 

D  WRITEC  6 , 8998  >  TIME 

8998  FORMAT<  '  TIME* ' ,  FI  0 . 5  > 

C 

C  GET  SIGMA,  V,  AND  TURBULENT  LENGTH  SCALE 
C 

100  CONTINUE 

CALL  SVLC SIGMA, V, EL) 

ALPHA  ■  V+DT/EL 

SIGMA  »  S1A1  .233*ALPHA**0.4976> 

D  WR I  TEC  6 , 8997  ) 

8997  FORMATC  '  SVL  '  ) 

C 

C  CALCULATE  POLES  AND  ZEROES 
C 

CALL  PANDZ<  PI  ,P2,P3,P4,Z1  ,,  Z2,Z3) 

D  WR I  TEC  6 ,8996  )  P 1 , P2 , P3 , P4 , Z 1 , Z2 , Z3 

8996  FORMAT<  '  POLES* ' , 4C 2X , D 1 5 . 9 >/ '  ZEROES*  ' , 3C 2X , D 1 5 . 9  )  > 

C 

C  CALCULATE  COEFFICIENTS 
C 

CALL  COEFFC  P 1 , P2 , P3 , P4 , Z 1 , Z2 , Z3 , A , B , C , D , F ,  XK 1  ) 

D  WRI TEC  6 , 8995  )  A, B, C, D, F, XK1 

8995  FORMATC  '  COEFF  A ,  B ,  C,  D,  F, XK1 * ', 6C 2X, E 1 5 . 9 > > 

C 

C  CALCULATE  RHS  DIFFERENCE  EQUATION  COEFFICIENTS 
C 

CALL  R I<  P 1 , P2 , P3 , P4 , A , B , C , D , F, R 0 , R 1 , R2 , R3 , R4  ) 

D  WRITEC 6, 8994 >  RO, R1 , R2, R3, R4 

8994  FORMATC'  RI= ' , 5C 2X, D1 5 . 9 > > 

C 

C  CALCULATE  DIFFERENCE  EQN  LEFT  SIDE  COEFFICIENTS 
C 

CALL  LIC PI ,P2,P3,P4, A, B, C , D , F , ELSO , ELS1 , ELS2 , ELS3 , ELS4  ) 
D  WRI TEC 6 , 8993  )  ELSO, ELS1 , ELS2, ELS3, ELS4 

8993  FORMATC'  LI*  ' , 5C  2X, D1 5 . 9  >  > 

C 

C  ITERATE 
C 
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YN  =  ELS1*Y3  -  ELS2*Y2  +  ELS3+Y1  -  ELS4*Y 0  +  R0*X4  -  R1*X3  +  R2*X2 
$,  -  R3*X1  +  R4*X0 

ITER  *  ITER  +  1 
SUM  =  SUM  +  YN 
SUMSQ  *  SUMSQ  +  YN*YH 
YBAR  =  SUM/ ITER 

Y2BAR  =  DSQRT<  SUMSQ/ 1 TER  -YBAR»*2> 

D  UR  I  TE<  6 ..  8900  >  Y3,  Y2,  Y1  ,  Y0,  X4 ,  X3,  X2 ,  X 1  ,  XO 

89  0  0  FORMAT <  '  Y3 , Y2 , Y 1 , Y 0= ' , 4< 2X , D 1 5 . 9  >/  ' , X4 , X3 , X2 , X 1 , X 0=  ' , 5<  2X , D 1 5 . 9  >  > 
YO  =  Y1 
Y 1  =  Y2 
Y2  =  Y3 
Y3  =  YN 
XO  =  XI 
XI  =  X2 
X2  =  X3 
X3  =  X4 

call  random  < sig, xb, xdp, xp, ix, iy  ) 

X4  =  ::dp 

D  UR  I  TEC  6 , 9998  )  X4 

9998  FORMAT<  '  RANDOM  NUMBER®  ' , F 1 0 . 6 > 

TIME  =  TIME  +  DT 
YP  =  YN 

URITE<6, 9999  >  T I  ME , YP 


c  write  data  to  plot  files 

c 

outv( 1  )=time 
outv<  2  >=yp 
do  1234  j=1 , icol 
nrec=nrec+1 

urite< 1 ,rec=nrec)  outv<j> 

1234  continue 

wri  teC3,  '<  e  1 3 . 7  )  '  >  yp 
irow=irow+1 
c 
c 

9999  FORMATS'  TIME® ' ,  E 1 2 , 6,  '  Y*',E15.9> 

IF< ITER  .LT.  NPTS  )  GO  TO  100 
UR  I  TEC  6, 9980)  NPTS , YBAR , Y2BAR 
9980  FORMATS  '  NPTS®", 18,'  YBAR®  ' , D 1 5 . 9 ,  '  Y2BAR®  ' , D 1 5 . 9 > 
c  IF<  LUPRT  . EQ .  8)  ENDFILE  8 

c 

wr i te< 2, '< 2 i5 > '  >  icol,irow 
wr  i  te<  2 , a8  >  '  )  'time' 
write*. 2, '< a8 > ' )  'turb' 
c 

STOP 

END 

c 

c 

SUBROUTINE  INITalC Y0, Y1 ,Y2,Y3,X0,X1 , X2 , X3 , X4 , T  IME , ITER) 
REAL*8  Y 0 , Y1 , Y2, Y3, X0, XI , X2 , X3 , X4 , t ime 
I NTEGER*4  ITER 
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X2  =  0 . DO 
X3  =  0. 0D0 
YO  =  0  .  DO 
Y 1  =  0 . DO 
Y2  =  O.DO 
Y3  =  O.ODO 
TIME  =  0. 
ITER  =  0 
RETURN 
END 


SUBROUTINE  SVL<  SIGMA,  V.,  EL  ) 
real*8  sigma, el, v 
V  =  100. 

EL  =  500. 

RETURN 

END 


SUBROUTINE  PANDZ< P 1 , P2, P3 , P4, Z 1 , 22, Z3 > 

COMMON/CFF/  SIGMA, V, ITER, AVK, EL, DT, PI 

REAL*8  PI ,P2,P3,P4, Z1 , Z2, 23 , C , P I , DT , P 1 P , P2P , P3P , P4P , 2 1 P , Z2P , Z3P 
real *8  avk , el , sigma, v 
INTEGER+4  ITER 

DATA  Z1 P/-1 . 46821 076D0/, Z2P/-5 . 3839321 43D0/, Z3P/- 0 . 6 1 2372436D 0/ 
8.  ,P1P/-1  .  02025059D0/ 

*,P2P/-1 . 6766  f  5? 1  DO/, P3P/-8 . 1  031 3370D0/, P4P/-1 .DO/ 

C  =  V/<  AVK+EL  1 
PI  -  P1P*C 
P2  =  P2P*C 
P3  *  P3P*C 
P4  *  P4P*C 
Z1  ■  Z1P*C 
22  =  Z2P*C 
Z3  =  Z3P+C 
RETURN 
END 


SUBROUTINE  COEFFCP1 ,P2,P3,P4,Z1 , Z2, Z3, A , B , C , D, F , XK 1  ) 
COMMON/CFF/  SIGMA, V, I  TER , AVK, EL , DT , P I 
REAL*8  P1,P2,P3,P4,Z1,Z2,Z3,A,B,C,D,F,XK1,DT,PI 
rea 1*8  avk,el,sigma,v 
I NTEGER*4  ITER 

XK 1  =  SIGMA*DSQRT<EL/< V*PI  )) 

XK1  =  XK1  *<  -*P  1  *P2*P3*P4  )/<  Z 1  *Z2*Z3  > 

A  -  XK1*<P1-Z1  >*<  PI -Z2  )*<  PI -Z3  )/<  PI *<  PI -P2  >*<  PI -P3  >*<  PI -P4  )  > 
B  =  XK1  *<  P2-Z1  >*<P2-Z2>*<P2-Z3)/<P2*<P2-P1  >*<  P2-P3  >*<  P2-P4  )  > 
C  *  XK 1 *<  P3-Z 1  >*<  P3-Z2  )*< P3-Z3  )/<  P3*C  P3-P1  >*< P3-P2 >*< P3-P4 > > 
D  *  XK1*<P4-Z1  >*< P4-Z2 >*< P4-Z3 >/< P4*< P4-P1  >*< P4-P2 >*< P4-P3  )  ) 
F  =  -XK1*Z1*Z2*Z3/<P1*P2*P3*P4> 

RETURN 
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SUBROUTINE  RI(P1 , P2,P3,P4,A,B,C,D,F,R0,R1 , R2, R3, R4  ) 

COMMON/CFF/  SIGMA, V, ITER, AVK, EL, DT, PI 

C0MM0N/EXPBL7  EX  1 , EX2, EX3 , EX4 , EX1 2, EX  1 3 , EX  1 4 , EX23 , EX24 
8,,EX34,EX123,EX124,EX134,EX234,X1234 
real*8  avk , e 1 , s i gma, v 

REAL+8  DT,PI  ,P1  ,P2,P3,P4,  A ,  B ,  C ,  [> ,  F ,  R 0 ,  R 1  ,  R2,  R3 ,  R4 ,  EX  1  ,  EX2 ,  EX3, 

&  EX4 ,  EX  1 2, EX1 3, EX  1 4 , EX23, EX24, EX34 , EX  1 23, EX1 24, EX1 34 , EX234 , X 1 234 
INTEGER+4  ITER 
EX1  -  DEXPCP1*DT> 

D  URITEC  6 , 9999  )  EX1 
9999  FORMAT<  '  EX1=',D15.9> 

EX2  =  DEXPC  P2  +  DT  ) 

D  URITEC  6 , 9998  )  EX2 

9993  FORMAT<  '  EX2=',D15,9) 

EX3  =  DEXP<  P3+DT  ) 

D  URITEC  6,9997  )  EX3 

9997  FORMAT<  '  EX3=',D15.9) 

EX4  =  DEXPC  P4*DT  > 

D  WRITER  6 , 9980  )  EX4 
9980  FORMAT<  '  EX4=  ' ,  D 1  5 , 9  !> 

EX  1 2=  DEXPC  CPI +P2 )*DT  > 

D  URITEC  6, 9996  )  EX12 

9996  FORMATC  '  EX1 2=  ' , D1 5 . 9  ) 

EX  1 3a  DEXPC  <  PI +P3  )*DT  ) 

D  URI TE<  6 , 9995  )  EX13 

9995  FORMATC '  EX1 3= ' , D1 5 . 9 > 

EX  1 4  =  DEXP<  C  P 1 +P4  )*DT  > 

D  URI TE<  6 , 9989 )  EX14 

9989  FORMATC  '  EX  1 4= ' , D 1 5 . 9 > 

EX23=  DEXPC  <  P2+P3  >*DT  ) 

D  URITE<  6, 9994  >  EX23 

9994  FGRMAT<  '  EX23=  ' , D 1 5 . 9 > 

EX24  -  DEXPC  <  P2+P4  >*DT  > 

D  URITE<  6, 9970  )  EX24 

9970  FORMAT<  '  EX24= ' , D 1 5 . 9  ) 

EX34  =  DEXP<  C  P3+P4  )*DT  ) 

D  URI TEC  6 , 9969  )  EX34 

9969  FORMAT<  '  EX34  =  ' ,  D 1  5  ■  9  > 

EX  1 23“  DEXPC  C  PI +P2+P3V+DT) 

D  URITE<  6 , 9993  >  EX123 
9993  FORMAT<  '  EX1 23= ' , D1 5 . 9 > 

EX124  ■  DEXPCCP1+P2+P4)*DT> 

D  URI TEC  6 , 9963  >  EX124 
9968  F ORMATC  '  EX  1 24= ' , D 1 5 . 9 > 

EX  1 34  ■  DEXPC  <  PI +P3+P4  >*DT  ) 

D  URITEC  6,9967  )  EX134 

9967  FORMATC  '  EX  1 34=  ' , D 1 5 . 9 > 

EX234  *  DEXPC  <  P2+P3+P4  )*DT  ) 

D  URI TEC  6 , 9966  >  EX234 
9966  FORMATC  '  EX234=  ' , D 1 5 . 9  > 

X 1 234  =  DEXPC  C  PI +P2+P3+P4  )*DT  ) 

D  URITEC  6, 9965  )  X1234 


9965  FORMATC'  X 1 234=  ' , D 1 5 , 9 > 
c 

R0»A  +  B  +  C  +  D  +  F  .-.j 

0  URITEC  6,9960)  RO  'f/A 
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9960 

c 

D 

9959 

c 


D 

9958 

c 


D 

995? 

c 


c 

D 

9992 


c 

c 


FORMAT (  '  R0=',D15,9> 

R 1  =  A+<  1  .  +  EX2  «■  EX3  +  EX4>  +  B*<  1  .  +  EX1  +  EX3  +  EX4  > 

&  +  C*< 1  .  +  EX  1  +  EX2  +  EX4  )  +  D*< 1  ,  +  EX1  +  EX2  +  EX3 ) 

&  +  F*<  EX1  +  EX2  +  EX3  +  EX4  > 

URITE<  6; 9959  >  Rt 
FORMAT< '  R1=',D15.9> 

R2  =  A*<  EX2  +  EX3  +  EX4  +  EX23  +  EX24  +  EX34) 

&  +  B*CEX1  +  EX 3  +  EX4  +  EX13  +  EX14  +  EX34) 

«.  +  C*<  EX  1  +  EX2  +  EX4  +  EX12  +  EX14  +  EX24) 

fc  +  D*<  EX  1  +  EX2  +  EX3  +  EX  12  +  EX  13  +  EX23) 

&  +  F*<  EX  1 2  +  EX  13  +  EX  14  +  EX23  +  EX24  +  EX34 ) 

URI TE<  6 , 9958  >  R2 
FORMAT<  '  R2= ' ,  D1 5 . 9  > 

R3  ■  A*<  EX23  +  EX24  +  EX34  +  EX234) 

&  +  B+CEX13  +  EX  1 4  +  EX34  +  EX134> 
fc  +  C*<  EX  1 2  +  EX  14  +  EX 2 4  +  EX124) 

&  ♦  D*<  EX  1 2  +  EX  1 3  +  EX23  +  EX123> 

&  +  F*( EX  1 23  +  EX  124  +  EX  134  +  EX234  ) 

WRITE<  6, 995?  ')  R3 
FORMAT <  '  R3=',B15.9> 

R4  =  A+EX234  +  B+EX134  +  C*EX124  +  D+EX123  +  F+X1234 
WRITER  6, 9992)  R 0, R 1 , R2 , R3 , R4 

FORMATC  '  RQ=',D15.9,  '  Rf»',t>l5.9,  '  R2=',D15.9,  '  R3= ' 

& ,  D 1  5 . 9 ,  '  R4  = ' ,  D 1  5 . 9  > 

RETURN 

END 


SUBROUTINE  LI< P 1 , P2 , P3 , P4 , A , B , C , D , F, ELSO , ELS1 , ELS2, ELS3 , ELS4 ) 
CONMON/'CFF/’  SIGMA,  V,  ITER  ,  AVK,  EL  ,  DT  ,  PI 

COMMON/EXPBL/  EX  1 , EX2, EX3, EX4 , EX  1 2 , EX1 3, EX1 4, EX23,  EX24 
&,EX34,EX123,EX124,EX134,EX234,X1234 
real*8  avk , e 1 , s i gma , v 

REAL*8  EX  1 , EX2 , EX3 , EX4 , EX  1 2 , EX  1 3 , EX  1 4 , EX23 , EX24 , EX34 , EX  1 23 
&  , EX  1 24 , EX1 34 , EX234 , X 1 234 , P I , DT 
REAL*8  PI ,P2,P3,P4,A,B,C,D,F,ELS0,ELS1 ,ELS2,ELS3,ELS4 
INTEGERM  ITER 
ELSO  =  1 .DO 

ELS1  =  EXJ  +  EX2  +  EX3  +  EX4 

ELS2  -  EX  1 2  +  EX  1 3  +  EX14  +  EX23  +  EX24  +  EX34 

ELS3  =  EX1 23*  +  EX124  +  EX134  +  EX234 

ELS4  =  X 1 234 

RETURN 

END 


BLOCK  DATA  EX 

COMMOH/CFFX  S I GMA , V , I  TER , AVK , EL , DT , P I 

COMMONS EXPBL/  EX1 , EX2 , EX3 , EX4 , EX  1 2 , EX  1 3 , EX  1 4 , EX23 , EX24 
&,EX34,EX123,EXJ24,EX134,EX234,X1234 
real*8  avk , el , s i gma, v 
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REAL*8  DT  ,  P I ,  EX  1  ,  EX2 ,  EX3 ,  EX4 ,  EX  1  2 ,  EX  1  3 ,  EX1  4 ,  EX23 ,  EX24 ,  EX34 ,  EX  1  23 
&  , EX  1 24 , EX1 34 , EX234 , XI 234 
INTEGER+4  ITER 
END 
c 
c 

subroutine  random  ( sig, xb , xdp , xp , ix, iy  ) 
rea 1+8  xp , a , xdp, sig, xb 
integer  iy, ix, j 
c 

a=G .  0 

do  30  j  =  1  ,  1  2 
c 

iy= ix*65539 
if < iy  )  10,20,20 
10  iy=iy+2147483647+1 

20  xp*iy 

xp=xp/2 1 47483647 . d0 


ix=iy 

30  a=a+xp 

xdp=< a-6. 0  >*sig+xb 
c 

return 
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